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ABSTRACT 

Experiments are carried out with various time and space 
differencing schemes applied to the barotropic primitive 
equations using both real data and a particular stream func- 
tion which is an analytic solution to the nondivergent baro- 
tropic vorticity equation. With both types of data, there 
were some significant differences in the forecasts produced 
by the various schemes. Replacement of the widely used leap 
frog (centered) scheme by others which eliminate or lessen 
some of its inherent errors at the expense of more computer 
time or storage appears to be justified at such time when 
computer capacity no longer restricts operational use of 
these more time consuming schemes. 
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I. INTRODUCTION 



In selecting a particular differencing scheme for appli- 
cation to numerical weather prediction, a major considera- 
tion is computer time and storage requirements. Therefore, 
the selected scheme is a compromise between verification ac- 
curacy and computer limitations. At present, the leap frog 
(centered) scheme is widely used in numerical weather predic- 
tion and will be used as a standard for purposes of compari- 
son. Present computer limitations and methods of verifica- 
tion tend to preclude seeking more accurate solutions based 
on mathematical considerations alone. The purpose of this 
study is to look beyond present computer capabilities and 
examine several other differencing schemes which reduce or 
eliminate some of the inaccuracies in the leap frog scheme. 
This scheme has three time levels and, therefore, a spurious 
computational mode arises from the finite differencing 
scheme. This problem can be eliminated by the application 
of two level time schemes such as the Euler backward or 
Runge Kutta methods which have on computational mode. A 
fourth order space differencing scheme is applied to reduce 
truncation error and improve phase speeds. Should any of 
the schemes produce a significant change in forecast rela^ 
tive to the leap frog method, it must be determined if the 
difference represents a sufficient increase in accuracy to 
warrant the added computer time and storage requirements. 
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II . PROGNOSTIC EQUATIONS 



The equations of motion and the continuity equation for 
a barotropic fluid in the momentum flux form are 



3(uh) 2 9 (uuh) 2 9 (uvh) ^ 

titt + m -;r~ ~ + m — ~ fvli = - mgh 



9t 



9x 



9y 



9(vh) . 2 9 (vuh) , 2 9 (vvh) , 

+ m -r — — + m -r — + fuh = - mgh 

9t 9x 2,, — 



9h , 2.9 (uh), 9 

•r— + m + — 

9 1 9x m 9y 



9y 
(vh) 



m 



9x 

9y 



m 



] = 0 



( 1 ) 

( 2 ) 

(3) 



where h is the depth of the fluid (or height of a pressure 
surface in the atmosphere), g is the force of gravity, f is 
the Coriolis parameter and u and v are horizontal velocity 
comp onen ts . 



A. SPACE DIFFERENCING 

1 . Second Order Scheme 

With u, V, and h given at each grid point, the 
differential equations may be approximated by 

(uh) . . _ 2 



uh , ^ . uh . . 

^ ij = -m'^/4d[(u. , .+ u..) ( ^ ^ 

9 1 < I \ ^ j ij ^ 'm m 



) 



(4) 



uh . - . uh . . 

- (u. , .+ u. .) ( ^3) 

1 - 1 , j ij m m 

+ (u. + u..)( ^3) 

i,j+l ij m m 

vh . , vh . , 

- (u + U )(— ^’3 + — 

1 » J “3- 1 j m m 



) ] + f (vh) 



i-J 



- mgh . . (h h. ^ .)/2d = F.. 

ij 1+1, J 1-1, J ij 

where j increases with x and i decreases with y, (See Fig, 
1) Similar expressions are used for the other two equa- 
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tions. This method of space differencing, devised by A. Ara- 
kawa (see, for example, II]), prevents nonlinear instability 
in the continuous time case (sometimes referred to as the 
semi-discrete equations) , 

2, Fourth Order Scheme 



account a greater number of grid points and reduces trunca- 
tion error approximates tendencies by 



B. TIME DIFFERENCING 

U represents u_h, or h; t is time; n is time step 

3 U 

number and F is in the following time differencing 

s chemes : 

1 . Leap Fro g 



Fourth order space differencing which takes into 







U 



n+1 



U 



n-1 



2AtF 



n 



2 . Runge Kutta 

AUi = At F[u", nAt] 

AU 2 = At F[u"+ AU^/2, (n+-|-)At] 
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4 
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At F[u" 


+ iUj . 


(n + 


1) 


At] 



u" + r[AU, + 2AU„ + 2AU. + AU, ] 
6 1 2 3 4 



3. Euler Backward 



* 
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u 



n 



At F 



n 



un+l_un 



* 

At F 



Denotes trial values. 
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III. BOUNDARY CONDITIONS 



A. A PARTICULAR STREAM FUNCTION 

The boundary conditions in this case are applied to an 
initial stream function which is an analytic solution to the 
nondivergent barotropic vorticity equation and after each 
time step during the forecast. 

1 . North-South Boundary Conditions 

For the 63x63 grid a wall exists between y=l and 
y=2 and also between y=63 and y=62 (See Fig. 1). No flux 
is permitted across this rigid boundary. This is accom- 
plished by setting v .= -v and v - = -v ^ thereby 

making the average 0 along the wall. The other para- 

meters are set equal across the wall. Summarizing, the 
boundary conditions are 

North South 



u 



4 ) 



X , 6 3 
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""x.l 


cvi 

X 
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X , 6 3 = 


u , ^ 


u , = 






X , 6 2 


X,1 


X , 2 


X , 6 3 


^x , 6 2 


-e- 

!-• 

II 


^x , 2 



where cf) is the geopotential. 

2 . East-West Boundary Conditions 

The East-West boundary conditions are cyclic as 
suggested by the function plotted in Fig. 1. Summarizing, 
the boundary conditions are 
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East 



West 



A.- ^ c. 'i 

l,y 62, y 63, y 2,y 

where A represents parameters v , u and d) 

x,y x,y ^x,y 

3 . Variation for Fourth Order Space Difference 

In order to apply fourth order space differencing 
and still maintain the boundary conditions of the original 
grid, it was convenient to create an additional outside row 
of points and relabel as a 65x65 grid, (See Fig. 2) The 
boundary conditions are summarized as follows: 
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B, ACTUAL 500 MILLIBAR DATA 

When the equations were integrated using real data, the 
^spongy” boundary conditions of the FNWC model were applied. 
After each time step, the prognostic results were modified 

by setting ^ - (i-k • ) uh + k ■ uh , 

vh = (l-k')A^ + k' vh „ 

h = (l-k')h + k' h- 
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where the 



subscript indicates initial values. Below lat- 



itude four degrees north set k’=4. Above 17 degrees north 
set k’=0. Between four degrees north and 17 degrees north, 
there is linear variation in k’ from zero to one. This re- 
stores the variables to their initial values south of four 
degrees north, while north of 17 degrees north, full varia- 
bility is permitted. There is no restoration north of 17 
degrees north. Between, there is partial restoration. 
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X*1 X'^3 




Illustrating the 63x63 square grid and 
boundaries as applied to a particular 
stream function. The dimensionless sym- 
bols X and y are used for notation con- 
venience . 



Figure 1 
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Illustrating the 65x65 expansion of the grid in 
Fig. 1 in order to apply fourth order space dif- 
ferencing over the same area. The dimensionless 
symbols x and y are used for notation convenience 

Figure 2 
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IV. GRID 



A. A PARTICULAR STREAM FUNCTION 

A 63x63 square grid as in Fig. 1 and as modified in 
Fig. 2 for fourth order space differencing is used to 
minimize the number of changes required to adapt programs 
to the FNWC 63x63 hemispheric grid. The boundaries in 
latitude extend from 30 to 90 degrees. Grid length is 
determined by dividing the geographic distance represented 
by this latitude difference into 62 equal parts. 



B. FNWC 63x63 HEMISPHERIC 

In conjunction with the actual 500 millibar data, the 
FNWC 63x63 hemispheric grid is used. This is a square re- 
presentation of a hemispheric polar stereographic grid in 
which the grid distance d is determined by 



d = 



m 



= dV( 



1 + sin 60 



0 



‘1 + sin Y 

where m is the map factor, d’ is the grid distance at 60^ 
latitude (381 kilometers), and sin y is defined by 

sin Y ■ ^^3.71 - „here =■ (I-I ) (J-J 

973.71 + r2 P P 



I,J = arbitrary grid point. 

I^,J^= grid point for the north pole = 32,32. 
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V, 



INITIAL CONDITIONS 



A. A PARTICULAR STREAM FUNCTION 

In order to obtain a smooth initial stream function 
field, an analytic solution to the nondivergent barotropic 
vorticity equation is used, namely 

I . .TT. dv 2 tT. 

^ sin ~ (y ~ “ (x-ct) 

y X 

where c is the phase speed and the wavelength in the x dir- 
ection and the half wavelength in the y direction are equal 
to 6ld (L^ = Ly = 6ld) • The initial stream function is 

, , . TT , dv 2 tt 

^ sin ~ (y “ j) — X . 

y . ^ 

Here, is determined by using the relation 






V L 
o X 

27T 



where v^ is the initial specified amplitude of the meri- 
dional wind speed derived from 

dip 2 TT , r . TT , dv 2 tT , 

♦o ~ 2' r~ ’'>■ 

X y X 

Using the initial field thus obtained, the linear balance 
equation 



v^(p = + Vf • Vip 



which in finite difference form may be written 
V^<})..= f . . V^ip. . + f V V f V iM E F.. 

IJ IJ A '• X X y y'*'-' xj 
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where 



’’x ■ ‘ ’i+l,j ■ * 'l-l.J 

The preceding equation may now be solved as a Poisson 

equation for the geopotential field. Thus, 

V^({) . . = F . . 
ij ij 

The initial guess for the relaxation procedure is 

where f is a mean coriolis parameter. 
Subsequent **guesses** are made according to the relation 



(}) = f \lJ 



(n+l)^ (n) ^ ^ (n) 

Y Y 



where 



- F. .) 

Ij Ij IJ 

and a is a relaxation coefficient. 

Using h = c()/g, computations for the initial height field 
are completed. Additionally, a constant height is added 
to each grid point value in order to eliminate negative 
heights, which would otherwise cause computational insta- 
bility. (See Fig. 3) 



B. ACTUAL 500 MILLIBAR DATA 



The initial height field is obtained from the FNWC 500 
millibar analysis. The linear balance equation is solved 
for the stream function as a Poisson equation using (j)/f 

as a first guess. 
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Differencing 
S cheme 



RMSE For Position 
(Nautical Miles) 



RMSE For Central 
Height (Meters) 



Leap Frog 


97 


30 


Leap Frog Fourth 


76 


31 


Order Space 






Runge Kutta 


93 


31 


Euler Backward 


88 
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Table of root mean square errors in position and central 
height of the three low centers at 24, 48, and 72 hours 
as forecast by various differencing schemes (on actual data) 
in comparison to the corresponding analysis. 



TABLE 1 




Initial height field as produced by a parti- 
cular analytic stream function on a 63x63 
grid. Only the central contours of the three 
features are illustrated. 

Figure 3 
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VI. RESULTS 



The results shown in the following figures were obtained 
from the computer subroutine mapping of the 63x63 height 
field obtained by the various differencing schemes using 10 
minute time steps. One grid squares are used in the figures 
simply as a convenience to provide a ready reference scale 
to compare positions of low centers. 

A. PARTICULAR STREAM FUNCTION 

Results for the forecasts made using a particular 
stream function for the initial values are summarized in 
Fig. 4 for central height and relative position. The fore- 
cast results for the individual differencing schemes are 
presented in greater detail in Fig. 8 through Fig. 12. 

B. ACTUAL 500 MILLIBAR DATA 

The results of forecasts made on three low systems from 
the 500 millibar FNWC analysis are summarized in Fig. 5, 

Fig. 6 and Fig. 7. Although only low systems are illustra- 
ted in this study, similar results were obtained from the 
movement and intensification of highs. The time frame of 
this study begins with the 1200Z March 24, 1966 500 milli- 
bar heights and terminates with the 1200Z March 27, 1966 
forecast. Low number one (See Fig. 5) was located over 
central Siberia. Low number two (See Fig. 6) and low num- 
ber three (See Fig. 7) were located over Southwestern 
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Canada and Labrador respectively. Root mean square error 
computations for position and central height as forecast 
by each differencing scheme when compared to the correspond- 
ing analysis are contained in Table 1. 

C, COMPUTER REQUIREMENTS 

Although peak program efficiency was not attempted dur- 
ing the application of the various numerical schemes, con- 
sistency in application was carefully maintained. Using 
the leap frog scheme as a standard, the leap frog fourth 
order space differencing and Euler backward schemes require 
approximately a 35 percent increase in run time. The Runge 
Kutta scheme requires approximately a 170 percent increase 
in run time and a ?5 percent increase in computer storage. 
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Figure 5 
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Forecast and analyzed relative positions and 
central heights for an actual data 500 milli- 
bar feature designated low number two. 

Figure 6 
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Forecast and analyzed relative positions and 
central heights for an actual data 500 milli- 
bar feature designated low number three. 

Figure 7 
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VII . CONCLUSIONS 



The various differencing schemes were applied to the 
smooth analytic solution in order that differences in fore- 
casts produced by these schemes would be more readily dis- 
tinguishable than with actual data. However, the choice 
as to resolution was probably too high in this study (61 
grid lengths per wavelength) and shorter wavelengths are 
yet to be tested. Differences produced in location and 
central values were not conclusive when compared to the 
results produced using actual data. A reduction in resolu- 
tion for the analytic initial field (approximately six grid 
lengths per wavelength) would possibly produce a more mean- 
ingful separation in the accuracy produced by the various 
schemes . 

The negligible increase in accuracy (four percent based 
on a root mean square error comparison) of the Runge Kutta 
scheme over the leap frog scheme does not justify the ex- 
tremely large increase in computer time and moderate in- 
crease in computer storage. Both the leap frog fourth 
order space and Euler backward differencing schemes show 
promise for improved accuracy with a moderate increase in 
computer run time (approximately 35 percent) . The leap 
frog fourth order space differencing scheme gave a 24 per- 
cent increase in position accuracy. The Euler backward 
♦scheme produced a nine percent increase in position accur- 



28 



acy and a 17 percent increase in central height. With both 
the Euler backward and the leap frog fourth order space 
schemes, the 35 percent increase in computer time could be 
reduced somewhat by more efficient programs. The small 
statistical sample of pressure systems examined in this 
study does not permit firm conclusions to be drawn regarding 
the various schemes tested. However, there is reason to 
believe that the fourth order space differencing of the 
nonlinear advection terms will reduce the phase error in the 
prediction of pressure systems and this method should be 
tested further. Williams [2] has shown that the pressure 
force need only be differenced with second order accuracy 
which permits a somewhat larger time step without computa- 
tional instability than when fourth order differences are 
used throughout. 
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One feature in an initial height field pro- 
duced from a particular stream function 
illustrating the central height and 390 
meter contour* 

Figure 8 
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Forecast height in meters of central values and 
positions of centers and 390 meter contours for 
a particular stream function by the leap frog 
scheme . 

Figure 9 
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Forecast height in meters of central values and 
positions of centers and 390 meter contours for 
a particular stream function by the Runge Kutta 
s cheme . 



Figure 10 
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Forecast height in meters of central values and 
positions of centers and 390 meter contours for 
a particular stream function by the Euler back- 
ward scheme . 



Figure 11 
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Forecast height in meters of central values and 
positions of centers and 390 meter contours for 
a particular stream function by the leap frog 
fourth order space differencing. 

Figure 12 
34 
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